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ABSTRACT 


An analytical framework was constructed to interchange ice resistance components from 
existing ice resistance calculation methods. Within this framework, the Lindqvist 
analytical method bending ice resistance component was substituted with a value 
obtained from a finite element static model of the ice sheet. A parametric study of this 
substitution was perfonned with error correction to within 12% over the entire range of 
parameter variation. The finite element ice sheet model was then damaged; and then the 
average ice bending resistance was obtained, substituted into the Lindqvist analytical 
method, and quantified as a change in the total ice resistance. Therefore, a hybrid ice 
resistance model was developed that accounts for the effect of damage to the bending 
resistance component, and enables further study of unconventional icebreaking methods. 
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I. INTRODUCTION 


A. ICEBREAKING HISTORY 

Prior to the 20th century, the first true icebreaker, Yermak, was commissioned 
and tasked with exploring the polar regions of the globe. The vessel had a displacement 
of 10,000 tons, a 10,000 hp propulsion system consisting of three screws aft with one 
screw forward, thickened hull plating, double hull construction, a ballast system to 
quickly affect roll/trim, and a rounded hull form [1], The summation of these design 
characteristics ensured this to be the first vessel to reliably transit multi-year sea ice, 
which it did from 1899 to 1963. The success of this ship ultimately defined the 
specifications for the IACS Polar Class of vessels, from which all modem icebreakers are 
derived [2], As a result of this tried and true formula for immensely powerful and 
specialized ice navigation platforms, the overall ship design and mechanical process to 
break ice is essentially the same today as it was back then. Despite this similitude, the 
icebreaker underwent several design iterations, and efficiency is considerably improved 
over 1899 as the concept was refined and matured. 

Substantial research was also conducted to develop systems to damage ice ahead 
of the icebreaker such that the ice failed with less effort from the ship. This concept was 
successfully used to aid early expeditions to the Antarctic, by fracturing the ice with 
explosives. Less violent means of damaging ice, such as water jet cutting mechanisms, 
were conceived to provide the benefit of the explosives with the practicality of a 
shipboard system. The Yermak served as the test bed for the first water jet cutters aboard 
icebreakers sometime after the idea was introduced in 1935 [3], The water jets at that 
time lacked the power density and peak pressures to provide any benefit to the 
icebreaker. At the time the Yermak was retired in 1963, follow on research attempted to 
physically model the effectiveness of an icebreaker with ice cutting systems [3]. The 
significant disparity between the failure characteristics of model and full scale ice 
resulted in criticism that the claims of this report were too optimistic in favor of ice 
cutting systems, and more specifically, water jet cutters [4], Nonetheless, the greatest 
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contribution of this work was the desire to use modeling methods to craft new strategies 
to damage the ice and allow the passage of a ship. 

From literature review, two themes following the 1963 study became apparent. 
First, the efficiency of the ice cutting systems received detailed study [4]-[8].Second, 
much research was done to adequately define the material properties of ice [9], [10], and 
apply this to reasonably model the icebreaking process [11], [12]. With minor exception 
[3], the research tended to ignore the application of these icebreaking systems in tenns of 
a grander icebreaking strategy. In other words, the specific damage pattern that causes ice 
to most readily fail in the presence of the icebreaker received little attention. Just as the 
summation of its parts cemented the Yermak’s title as the first true icebreaker, the 
effectiveness of the ice breaking / ice cutting system relies on the synergy between the 
available tools and their artful manipulation. 

B. STUDY OBJECTIVE 

This study combines the advantages of existing ice resistance calculation methods 
into a single hybrid method, such that the steady state velocity of an icebreaker is 
correlated to the available thrust; and include cases in which the ice sheet possesses some 
artificial damage. In addition, this study intends to demonstrate the utility of this hybrid 
model by rationalizing ice resistance estimates against other modeling techniques, full 
power test results, and realistic performance of modern ice cutters. 
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II. BACKGROUND 


A. DEFINITIONS 

To more appropriately convey the subject material, all pertinent variables are 
hereby defined in Tables 1-5. All equations obtained from literature or developed are 
based on these definitions to aid consistency in presentation. Table 1 provides the unit 
vectors for two body-fixed coordinate frames. The icebreaker reference frame is given in 
Cartesian coordinates. The icebreaking surface frame is defined in two dimensions and is 
coupled to the icebreaker reference frame. The parameters and dynamics for the 
icebreaker are provided in Table 2. All parameters are typically known after the initial 
design stage. The icebreaker body has six degrees of freedom; three translational 
motions, and three rotational motions. This study emphasizes forward translation of the 
icebreaker only within this context, although follow on study should include additional 
motions. The parameters of sea ice including dimensional, material property, and material 
failure considerations are provided in Table 3. All physical parameters and variables for 
existing analytical ice resistance methods are included in Table 4. All parameters 
affecting the development and error of the hybrid method for average ice bending 
resistance are included in Table 5. An illustration of the reference frames, hull 
parameters, and motions are provided in Figure 1. 
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Table 1. Body-fixed coordinate system conventions. 


Reference Frame 

Coordinate System 

Unit Vector 

Definition 

Icebreaker Body-Fixed 

X 

longitudinal 


Y 

transverse 


Z 

vertical 

Hull-Ice Contact Surface 

n 

normal 


t 

tangential 


Table 2. Icebreaker parameters and motions variable definitions. 


Icebreaker 

Parameters 

Variable 

Name 


A 

displacement 


M 

mass 


Ps 

shaft power 

Dimensions 

L 

length at waterline 


B 

breadth 


T 

draft 

Bow Angles 

a 

waterline entrance 


3 

flare 


Y 

stem 


? 

effective ice contact 

Dynamics 

Variable 

Motion Affected 

Translation 

X 

surge 


y 

sway 


z 

heave 

Rotation 

0 

roll 


4> 

pitch 


43 

yaw 


4 
































Table 3. Sea ice parameters. 


Sea Ice 

Parameters 

Variable 

Name 

Dimensions 

d 

grain diameter 


H 

thickness 


Ic 

characteristic length 


1 

cusp length 


w 

cusp width 


R 

cusp radius 


Ro 

cusp radius reference 


n 

cusp radius angular position 

Properties 

E 

Young's modulus 


pice 

density 


V 

Poisson's coefficient 


P 

friction coefficient 


va 

air concentration 


vb 

brine concentration 


Q 

temperature 

Material Failure 

ob 

bending stress 


ob,o 

uncorrected bending stress 


oc 

crushing stress 


rs 

shear stress 


ol 

principle stress 


o2 

principle stress 


do/dt 

stress rate 


Sol/Sx2 

stress gradient 


C 

stress rate constant 


D 

stress gradient constant 


K 

stress gradient limit 
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Table 4. Physical parameters and resistance component variable definitions. 


Resistance Calculation 

Parameters 

Variable 

Name 


pw 

sea water density 


Fp 

propulsion force 


Fv 

vertical ice loadingforce 


g 

gravitational acceleration 


n 

number of cusps 


V 

velocity 


Vi 

initial velocity iterate 


Vi+1 

product velocity iterate 

Resistance 

Rb 

bending resistance 


Rc 

stem crushing resistance 


Rice 

ice resistance 


Row 

open water resistance 


Rr 

rotation resistance 


Rs 

submergence resistance 


Rsf 

skin friction resistance 


Rw 

wavemaking resistance 


Table 5. Hybrid method parameters for average ice bending resistance. 


Finite Element Bending Resistance 

Parameters 

Variable 

Name 


A 

FEM proportionality constant 


b 

error slope-intercept 


F 

distributed ice loadingforce 


m 

error slope 


N 

parameter power constant 


P 

substitute parameter 
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B. LITERATURE REVIEW 

1. Sea Ice Properties 

Accurate ice resistance models require knowledge of sea ice material properties 
and failure mechanisms in the presence of an icebreaking ship. However, the material 
properties of sea ice are difficult to assign because many studies identify variations in ice 
samples based on age, chemical composition, temperature, scale, and testing method [9]. 
For example, in the Lindqvist 1989 study [13], nearly one order of magnitude variation in 
the flexural failure stress of ice was observed during in-situ testing. Therefore, the 
material properties obtained from literature are approximate quantities, and likely require 
refinement based on further specification of these variables. 

Sea ice is fonned from accumulated snow at the surface of seawater, and has a 
microstructure that differs considerably from freshwater ice. Sea ice typically assumes a 
hexagonal column structure with anisotropic material properties. Unlike freshwater ice, 
which forms a uniform crystalline structure, the grains in sea ice orient in several 
directions, lending to the formation of an ice sheet with varying thic kn ess. Additionally, 
brine channels form between the grains, forming a heterogeneous solid solution within 
the ice sheet, and giving the ice viscoelastic properties. Repeated tests indicate either 
ductile or brittle material failure depending on the ice temperature, composition, and 
stress rate [9]. 
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a. Density 


•5 

The density of sea ice is roughly approximated as 910 kg/nr [9]. 

b. Young’s Modulus 

Given the viscoelastic behavior of sea ice, Young’s Modulus varies from 
0.3 GPa-10 GPa for static loads, to 6 GPa-10 GPa for dynamic loads [9]. 

c. Poisson’s Ratio 

The Poisson’s ratio is approximately 0.295 [9], although it is also 
presented as 0.3 [9]. 

d. Flexural Strength 

(1) Temperature 

The flexural strength is related to the temperature of the ice for 
temperatures -0.4 C > Q > -18.5 C [9]. 

cj b =4.7-0.960-0.3 IQ 2 (2.1) 

(2) Composition 

The expression for flexural strength as a function of composition is 
shown below, where v a is the concentration of air, and Vb is the brine concentration [9]. 

°b = °b,o( V a’ V b) ( 2 - 2 ) 

(3) Stress Rate 

The flexural strength of the ice has no dependence on stress rate in 
the range 49.0 kPa/s < < 981 kPa/s. For flexural stress rates above 981 kPa/s, the flexural 
strength is obtained as a function of the stress rate [9]. 

. da _ 

v b =v h ,o + c{n — ( 2 - 3 ) 

dt 
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(4) Testing Method 

Depending on either a small beam or ring tensile test for ice, the 
tensile stress of ice was measured at 0.785 MPa and 1.77 MPa, respectively [9]. 

(5) Grain Diameter 

The flexural stress of ice is approximated as a function of the grain 
diameter, typically between 1 cm and 2 cm [9]. 


Where the grain diameter, d is given in centimeters, and the flexural stress, Ob is in MPa. 

e. Shear Strength 


The shear test estimates the shear failure stress at 2.06 MPa [9], 


f. Compressive Strength 

Compressive failure of the ice varies by orientation of the ice crystal, 
having a maximum of 8.34 MPa for cores loaded on the principle axis, and a minimum of 
3.43 MPa for cores loaded perpendicular to this principle axis [9]. 


g. Stress Gradient 


In addition to these three failure mechanisms, sea ice fails at a stress 
gradient threshold [10]. 


K 2 =(a l -a 2 ) 2 +<J, 2 +<J 2 2 ~D 


' da ^ 2 

\dx 2 j 


(2.5) 


This equation modifies the von-Mises plane stress failure criterion to 
include the stress gradient term, the partial derivative of the principle stress with respect 
to the perpendicular axis. 


2. Icebreaking Process 

Icebreakers induce failure of the ice by three observed mechanisms; crushing, 
shearing, and flexure. In most icebreaking operations, all three of these mechanisms are 
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present to some degree. Crushing and shear failure are localized at the hull-ice contact 
area, whereas flexural failure occurs at some distance from the icebreaker. The general 
pattern of damage generated by an icebreaker is a crushing zone at the stem of the 
icebreaker, with the formation of characteristic semi-elliptical cracks, kn own as cusps, 
from the centerline to the extreme breadth, caused by flexural failure. The cusps may also 
exhibit secondary crack formation from the contact surface radially along the length of 
the cusp [13]-[16]. 


Cusp Formation 



As a result of the material properties of sea ice, this pattern of cusp formation is 
dependent on icebreaker dynamics, as observed in laboratory and scale-model testing 
[15]. In certain conditions, the ice cusps occur predictably at specific locations along the 
bow. However, past research cannot conclude whether icebreaker dynamics are beneficial 
or detrimental to the overall ice resistance [15]. In addition, the interplay of failure 
mechanisms is dependent on the dimensions of the ice due to stress gradient limitations 
[10]. Upon failure of the ice, the icebreaker clears the ice underneath and to either side. In 
this stage, the icebreaker both rotates and submerges the broken pieces, providing some 
resistance to the ship. However, as the broken pieces are no longer providing structural 
resistance from the ice sheet, the thrust exceeds the resistance forces, and the icebreaker 
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regains momentum. To maintain transit speed through the ice, all momentum lost by the 
loading of the ice must be restored within a distance equal to the length of the broken 
pieces in the direction of motion, before the icebreaker once again encounters the ice 
sheet [16]. 


3. Model Techniques 

Since icebreaking ships traditionally represent significant capital investment, 
modeling is conducted to allow iterative performance improvements to the intended 
design. Three types of modeling were developed, each with inherent advantages and 
disadvantages depending on the application. The three types of models are; scale models, 
analytical models, and numerical models. 

a. Scale Models 

Similar to tow ta nk testing, this modeling technique incorporates a layer of model 
ice floating on top of the water. The desired hull form is constructed, instrumented, and 
physically tested while breaking model ice. This technique is advantageous because 
resistance data is obtained directly, according to scale, and involves physical failure of 
the sea ice. The major caveat to this advantage is that the stress gradient, stress rate, and 
complex composition of sea ice disrupt the accuracy of the scale results. To more 
accurately simulate the full scale icebreaking process, the model ice is either “doped” by 
adding chemicals to affect the material properties, or is manufactured similarly to natural 
sea ice by allowing artificial snow to freeze on the tow ta nk surface [15]. Despite realistic 
interpretation of full scale icebreaking operations, the disadvantage to this method is that 
it is not practical at the initial design stages because of cost, time, and facilities 
availability. 


b. Analytical Models 

In contrast to scale modeling, analytical models [13], [14] are useful to 
calculate ice resistance at the initial design stages because of their simplicity. For ships 
following the general specifications of a modem icebreaker, existing analytical models 
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are the most practical means to determine ice resistance. Since analytical models are 
constrained by a certain set of assumptions for a particular method, a major disadvantage 
of a semi-empirical analytical model is that results are scaled with respect to existing 
similar ships analysis, and may not provide reasonable performance estimation for more 
radical ship designs. Even with a standard hull form, in the case of the ice cutting ship, 
current analytical methods are constrained to local ice failure, and cannot account for the 
effect of damage to the ice sheet at some distance from the icebreaker. However, one trait 
of analytical methodology is the decoupling ice resistance into various constituents. For 
example; open water resistance, bending resistance, crushing resistance, rotational inertia 
resistance, and submergence resistance may be solved for separately, according to 
approximate system parameters, and scaled to match empirical data. Therefore, it is 
feasible that some modification of particular resistance constituents within the analytical 
method may allow inclusion of ice sheet damage to the overall ice resistance calculation. 

c. Numerical Models 

While finite element models have the potential to replicate the physical 
phenomena of breaking ice, these models encompass a wide range of programs with 
varying complexity and accuracy. Numerical methods have already implemented variable 
parameters and dynamics affecting icebreaking processes [17]-[19]. However, the 
inherent imperfections and nonlinearity of the solution of the numerical model requires 
validation with empirical data across the entire range of expected operations, adding to 
the cost of the technique. Although numerical models may provide insight into new 
radical icebreaker designs at marginal expense, they may also contain significant 
unforeseen errors that must be accounted for by other means. 

4. Cusp Geometry 

Several research efforts attempt to relate ice cusp dimensions to the ice 
loading because of its significant influence on ice resistance predictions. As a result, three 
types of ice cusp approximations were obtained from literature, with varying application 
and accuracy describing this form of damage realistically. In order of complexity, ice 
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cusps were analyzed according to one-dimensional, two-dimensional orthogonal, and 
two-dimension radial coordinate systems as shown in Figure 3. Within this spectrum of 
analysis, the ice is either uniformly loaded along the edge, or experiences point loading. 
The relative accuracy of each depends on the geometry of the hull. Since uniform loading 
is expected of bows with flat contact surfaces, such as the Thyssen-Waas bow, the 
uniform loading analysis can reasonably describe the icebreaking pattern. For other bow 
forms with some curvature, however, point loads are a more realistic assumption. In all of 
these cases, the ice is assumed as a homogeneous and isotropic sheet of uniform 
thickness with perfectly vertical edges on a uniform elastic foundation. 



a. One Dimensional 

The least complex form of ice cusp analysis assumes the ice sheet as a 
semi-infinite plate loaded along a single free edge. In this form, the ice is assumed to fail 
at a set distance from the free edge of the plate, proportional to the characteristic length 
of the plate, Ic [13]. 
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This distance from the edge to the location of failure is referred to as the cusp length. Since 
this method fails to account for the cusp width, two issues arise that affect its application. 
First, it is difficult to derive a realistic physical association between the assumed distributed 
load required for cusp failure and the point load it represents. Second, the cusp length 
detennined by this method is likewise restricted to one-dimensional applications. 
Therefore, such analysis is best suited to cases exhibiting: 1) distributed loading of the ice 
sheet, and 2) the width of the distributed loading is much greater than the cusp length. For 
the majority of icebreakers, however, this method is not reliable. 


b. Two Dimensional Orthogonal 

There are two variations of this method, constituting the progression from a 
distributed load to a point load assumption. This is accomplished by assuming specific 
failure shapes according to prescribed geometries, in which the orthogonal cusp width is 
related to the cusp length. The cusp shapes are defined by either rectangular or triangular 
elements, as shown in the figure. In all cases, the cusp length is determined by the same 
method as the one-dimensional analysis. It is important to note that several analytical 
icebreaking resistance methods combine the rectangular and triangular elements, 
therefore, accepting a significant departure from both the natural icebreaking process and 
analytical description. 

(1) Rectangular Element 

This variation is similar to the one-dimensional method, except that the cusp 
width is related to the cusp length. This allows the frequency of ice cusp formation to be 
determined, thus improving the accuracy of the total ice load approximation, with 
marginal increase in calculation. However, as in the one-dimensional case, it is difficult 
to assign an appropriate loading force that bears physical compatibility. 

(2) Triangular Element 

This variation relates the length of the cusp as sides of an isosceles triangle to an 
arbitrarily determined base, the cusp width. This allows the cusp to be solved in the same 
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manner as the rectangular cusp, with exception of variable cross-section along the cusp 
length. However, since the free edge is infinitely small, the distributed force acting on 
this representative cusp may be approximated as a point load. In addition, these triangular 
cusps may be joined side by side to form a cusp resembling a semi-circle experiencing a 
point load at the center. If an infinite number of triangular elements are used, the error 
between the cusp shape and triangular representation can be effectively eliminated. 
Therefore, this variation allows both versatile and somewhat realistic description of the 
cusp. However, the solution complexity limits it to numerical solvers [20], 

c. Two Dimensional Radial 

The most complex, and likewise, accurate methods for determining the 
relationship between the ice loading force and semi-elliptical cusp dimensions 
appropriately analyze the system in a two-dimensional radial system. Similarly to the 
triangular cusps, the ice is loaded at a single point. Since the relationship obtained from 
literature has empirical basis, it is much more easily calculated upon detennining the cusp 
variables. The caveat to this, however, is that a numerical solver is required to determine 
all variables [18]. 
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III. LINDQVIST METHOD 


A. REASON FOR COMPARISON 

The Lindqvist method for calculation of ice resistance is a semi-empirical 
analytical model that is referenced by several modern ice resistance methods as a 
standard for comparison [16], [21], [22]. The Lindqvist method likewise serves as a 
benchmark for this particular study, although follow on research indicates that the 
method likely has some conservative error. Therefore, the hybrid ice resistance method is 
developed out of respect to the Lindqvist method, adopts components of this method to 
reduce error in comparison, and attempts to make improvements to experimental error. 

B. ASSUMPTIONS 

Several assumptions were made in the development of the Lindqvist method that 
affects the accuracy of the solution. These assumptions fall into two general categories, 
those that are in place to eliminate nonlinearity from the method, and those that provide 
rough approximation for poorly understood quantities. The empirical basis and linearity 
of the method allow it to be reasonably accurate, despite room for improvement within 
the individual assumptions. The assumptions of the Lindqvist Method are as follows: 

1. Linear Assumptions 

• Ice resistance increases linearly with velocity. 

• The vessel is motionless with exception of constant forward translational 
velocity. 

• The hull is approximated as a plane, developed from average hull angles, 
exerting normal and tangential force components on the ice at the 
waterline. 

• Upon contact, the ice fails instantaneously. 

• The failed ice pieces do not provide rotational inertia resistance back to 
the hull, nor are they affected by the viscosity of the water. 

• The ice does not exhibit viscoelastic, stress rate, or stress gradient 
properties. 
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2. Quantity Assumptions 

• Broken ice provides frictional resistance evenly distributed along 70% of 
the underside of the hull once submerged. 

• The resistance force is independent of the crushing failure stress. 

• The shear failure stress equals the bending failure stress. 

• Young’s modulus for ice is 2 GPa. 

• Poisson’s ratio for ice is 0.3. 

• Friction coefficient during hull-ice interaction is 0.1 for hulls with low- 
friction coatings and 0.16 otherwise. 

• Both the length and width of a failed ice cusp are equal to one third of the 
characteristic length of an edge-loaded semi-infinite plate with uniform 
elastic support. 


C. RESISTANCE COMPONENTS 
1. Stem Crushing 


The equation derived to approximate the stem crushing resistance is provided 
below. The equation is multiplied by an approximate proportional constant obtained by 
instrumenting the bow of the icebreaker with strain gauges, and taking transient 
measurements with reference to the ice bending failure stress and thickness. It is 
important to note that this resistance is not influenced by the dimensions of the 
icebreaker. 


Rc = 0-5 


<Ju 
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cos Y 

tan y + n -- 

cos Q 


1 -// 


siny 

COStT 


(3.1) 


where the effective ice contact angle between the hull and the ice, is defined as a 
function of the waterline entrance and stem angles. 


C, = tan 1 


tan y 
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(3.2) 
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2. Bending 

The bending resistance equation assumes the failure of a non-integer number of 
ice cusps with dimensions proportional to the characteristic length of the ice sheet, 
modeled as a semi-infinite plate supported on an elastic foundation. 

J n cos y (. 1 

tan^ —- - — 1 +- 

v sin a cos cos^ 

3. Submergence 

The submergence resistance component provides the frictional contribution of the 
broken ice flows as a function of the approximated underwater surface area of the hull 
and submerged depth. This formulation incorporates the snow thickness resting on top of 
the ice because this is also submerged. This equation does not explicitly consider the 
influence of the added mass of water surrounding the ice chunks as they are submerged. 

K=(P,-p,„)g(H,«+H„)B (3.4) 

yn + 21 J 

where 
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D. EMPIRICAL VELOCITY RELATION 


As Lindqvist observed, measured ice resistance appears to increase linearly with 
velocity from some nonzero resistance value as expressed below. 


K,=(K + R t ) i + i-4 


+ ft 1+9.4- 


The Lindqvist empirical velocity relation accounts for this observation by setting the 
zero-velocity resistance as the sum of the three resistance components, and increasing the 
value of each component linearly according to velocity, and scaling these values 
according to non-dimensional constant values. Thus, this velocity relation implies that 
two criteria must be met. First, assuming a linear relation between ice resistance and 
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velocity as observed, the three resistance components must equal the predicted zero- 
velocity slope-intercept of the empirical resistance data. Second, the contribution of each 
resistance component must be scaled appropriately with velocity, such that the sum of 
three resistance components equals the total resistance for a range of vessels. 

E. INTRINSIC VALIDATION 

The Lindqvist method was validated against both the calculated and empirical 
resistance data extrapolated from the original study. This comparison was necessary to 
ensure that the method was implemented properly into this study, since the Lindqvist 
resistance is typically presented as the total resistance, without specifying the 
contributions of the underlying resistance components. As desired, the resistance values 
calculated by Lindqvist were successfully reproduced for each vessel in the original 
study. Therefore, the individual resistance components analyzed by Lindqvist during the 
formation of his method were likewise reproduced, as shown in Figure 4. 



Figure 4. Relative influence of Lindqvist resistance components on the icebreaker 

Vladivostok with respect to velocity. 
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This analysis provided several insights into limitations of the method. Although 
the final empirical velocity relation achieves reasonably accurate calculation of the total 
ice resistance, there may be substantial errors in the individual resistance components. 
For example, the average resistance component values were obtained and compared with 
the total resistance versus vessel velocity. This indicated that although ice resistance 
appears to increase linearly with velocity, the contributions of the stem crushing and 
bending resistance components increase, while the submergence resistance component 
decreases. In addition, the relative contributions of each of these three resistance 
components were compared. Averaging from 0 m/s to 5 m/s, ice submergence provides 
approximately 40% to the total resistance, stem crushing provides 35%, and bending 
contributes 25%. Therefore, these percentages allow initial approximations of the ice 
resistance for icebreakers designed to eliminate any one component. 



Figure 5. Extrapolated raw data for the icebreaker Vladivostok comparing the original 
Lindqvist ice resistance calculation with a calculation using the method and 

the best fit line of the data. 
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Using the extrapolated data, and assuming a linear resistance velocity 
relationship, a best fit line was drawn from the raw data, as shown in Figure 5, since this 
was not provided in the original study. Overlaying the Lindqvist calculation over the raw 
data and best lit line revealed that although the Lindqvist data reasonably approximates 
the resistance, the zero-velocity resistance, nor the slope of the resistance line are 
properly matched by the method. This was expected from the assumptions of the method 
and analysis over a range of vessels, however, these errors may also stem from the 
absence of a pertinent resistance component, such as rotation of the ice. 

F. SENSITIVITY ANALYSIS 

Successful implementation of the Lindqvist method enabled sensitivity analysis to 
determine the most influential variables on total ice resistance with regard to the 
individual resistance components, located in Appendix A. Table 5 lists the parameters 
used for the sensitivity study of each ice resistance component. 

Table 6. Resistance component influential parameters 
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Lindqvist acknowledges that some portions of the method require improvement, 

therefore this analysis contains errors. Since Lindqvist also indicates that his method is 

more accurate for larger vessels, the largest vessel of the original study, Vladivostok, was 

selected to mitigate these errors. As expected, the sensitivity analysis determined that the 

ice resistance is highly dependent on the vessel gross dimensions, bow angles, and the ice 
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dimensions. However, the Lindqvist method proved highly sensitive to the bending stress 
of the ice and the density differential between sea water and sea ice. This dependence on 
environmental variables requires that ice resistance calculations adequately describe the 
intended operating environment. Although the Lindqvist study described the bending 
stress for all cases, it did not describe the density differential. After several iterations, the 
best error obtained from adjusting the density differential was to within 10% of the 
original study. Therefore, the exact value of the density differential was not reproduced, 
but was close enough to proceed with analysis. 
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IV. HYBRID METHOD DEVELOPMENT 


The proposed hybrid method attempts to evaluate both conventional and other 
icebreaking strategies by achieving synergy among all modeling methods. The creation of 
this tool is motivated by the desire to perform objective icebreaking resistance analysis 
under a common framework, as shown in Figure 6. This framework is intended to 
provide the context to combine resistance components from multiple sources to 
determine the overall resistance. Comprehensive analysis of this method is not possible 
within the scope of this thesis because of its many variations across the range of all 
resistance approximation methods. Therefore, the proposed system requires further 
analysis, refinement of individual resistance components, and their relations. 

A. PROPOSED FRAMEWORK 

1. Resistance Components 

Research was conducted to detennine all analytical resistance components from 
literature, listed below. 

• Bending 

• Rotation 

• Skin Friction 

• Stem Crushing 

• Submergence 

• Wave Making 

2. Icebreaking Functions 

Separation of resistance components was accomplished by first subdividing the 
icebreaking process into three necessary functions; breaking the ice, clearing the ice, and 
making headway. 


a. Icebreaking 

The stem crushing and bending resistance components are included within 


the icebreaking function because these resistance components are only applied as the 
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icebreaker comes into contact with the ice sheet although these resistance components are 
often summed as steady values, they are actually cyclical depending on the failure of the 
ice due to the failure of the ice extending some distance ahead of the icebreaker. The 
stem crushing component provides more steady resistance to the icebreaker since the 
stem is the first part of the icebreaker to contact the ice, and is relatively unaffected by 
cusp failure. The bending component cyclical frequency, on the other hand, is highly 
dependent on the dimensions of the ice cusps. 

b. Ice Clearing 

The clearing of broken ice involves rotation and translation of the broken 
chu nk s. This function is accomplished by the rotation and submergence resistance 
components existing in literature. The rotation ice resistance component is a function of 
both the rotational moment of inertia of the failed ice cusp and the icebreaker velocity. 
The submergence resistance component incorporates both the translation of the broken 
ice under the water and the hull-ice contact friction, however, these may also be 
separated. 


c. Hydrodynamic Resistance 

Progressing relative to the ice sheet also involves making headway relative 
to the water beneath the ice. This relative motion creates drag on the icebreaker body in 
the form of frictional and pressure drag. The skin friction resistance component is directly 
proportional to the surface area of the hull in contact with the water. The pressure 
differential generates waves through this medium at the water-air interface. This drag 
component is defined in literature as the wave making resistance. Since the skin friction 
and wave making resistance components are traditionally of great importance ship 
design, they are typically obtained as a coupled resistance from scale model tow tank 
testing, but may be accounted for separately by a wide range of other methods. 
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Figure 6. Proposed icebreaker resistance framework. The corners of the triangle 

represent the three icebreaking functions. The Venn-diagram illustrates the 
icebreaking resistance components obtained from literature. The arrows 
indicate potential relationships between resistance components. 


B. ANALYTICAL VALIDITY 
1. Resistance Coupling 

The success of this method is governed by the validity of the assumption that the 
icebreaking process is comprised of discrete resistance components. This assumption is 
not particularly radical since the resistance components within the framework were all 
obtained from literature, and have been implanted into existing resistance models 
successfully [13], [16], [18]. Flowever, Figure 6 also identifies cases, both described in 
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literature and newly added, in which some resistance elements are coupled. In the broader 
sense, this illustrates the interdependence of the three icebreaking functions. The first 
such coupling was observed in previous research, that hydrodynamic forces resulting 
from the motion of the icebreaker contributed to failure of the ice by affecting the 
buoyant force supporting the ice sheet [16]. This has also been observed during high 
speed icebreaking by hovercraft [6]. The remaining three instances of coupling are 
speculated in this study. The first speculates that the bending and rotation resistances are 
linked by the geometry of the failed cusp. The second potential relationship links the 
submergence and hydrodynamic skin friction components by the finite wetted hull area. 
The third assumes that there is some connection between stem crushing and bending 
components based on the dimensions and form of the bow. 

2. Coupling Mitigation 

This lack of infonnation provided the motivation to speculate which parameters 
influence this coupling such that analytical errors are reduced or expected. 

a. Velocity 

The wave making and skin friction resistance values become trivial as 
velocity approaches zero. Since these values also contribute to coupling of the bending 
and submergence components, respectively, it is reasonably assumed that analytical error 
is reduced as velocity approaches zero. The zero-velocity assumption does not have 
serious implications on icebreaker operations since it produces the operating limits. 

b. Hull Geometry 

The stem crushing, bending, and rotation resistance components are 
effectively negated as the hull-ice contact angle approaches zero. However, these values 
remain coupled. Since the only design parameter not shared by all three of these 
components is the vessel breadth, this value may be adjusted to reduce expected 
analytical error. The stem crushing resistance component is independent of vessel 
breadth. Therefore, the stem crushing is effectively de-coupled from the bending 
resistance in cases where the vessel is narrow, or has an extremely large breadth. Two 
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vessels that may provide further insight into this assumption at each extreme are the E- 
Craft and the oblique icebreaker concept, respectively. 

C. RESISTANCE COMPONENT SUPERPOSITION 

Despite the validity of assuming truly decoupled resistance components, the 
models presented in literature approximate icebreaking resistance as the sum of either 
average or discrete values over the icebreaking cycle. If a steady state resistance value is 
desired at a specified velocity, the resistance components obtained from either means 
may be used interchangeably according to the work-energy principle. 



In addition, these values may be interchanged to obtain a rough resistance estimate for 
nonlinear systems. 

1. Conceptual Test 

To test this concept, the Lindqvist calculation was performed for the icebreaker, 
Vladivostok. In the first test, it was assumed that the average value of all resistance 
components were present through the entire duration of the icebreaking cycle. The 
steady-state velocity of the icebreaker was obtained via application of the work-energy 
principle. In the second test, the bending resistance component was assumed to occur 
only upon contact with the ice sheet until failure of the ice cusp, with all cusps failing 
simultaneously. To further test this concept, the second test was also modified such that 
that individual cusps were assumed to fail at equal intervals. 

a. Averaged Ice Resistance 

Upon applying the work-energy method presented in Equation 4.1 to the 
transient numerical model, the icebreaker achieved steady-state velocity as displayed in 
Figure 7. This method includes the mass of the icebreaker, and therefore, more accurately 
represents the velocity profile of an icebreaker. In this averaged ice resistance case, all 
resistance values are assumed to be present at every velocity iterate. In addition, the ice 
resistance increases proportionally with velocity. Since the hydrodynamic resistance was not 
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included in this model, the increase in ice resistance causes the convergence of the velocity. 
In this particular case, the steady state velocity converged to approximately 0.325 m/s. 


Vladivostok: Continuous Cusp Failure Velocity 
Profile at R ice = 600 kN 

0.35 ; 



Figure 7. Continuous cusp failure transient velocity profde. 


b. Cyclical Bending Resistance 

The steady-state velocity profde for simultaneous and single cusp failure 
was obtained using Equation 4.1, appearing as a saw-tooth pattern as expected from 
literature [13]. Figure 8 illustrates the contrast between assuming simultaneous and 
individual cusp failure. This model differs from the average ice resistance model because 
it has the added dependence on the absolute distance traveled, along with velocity 
dependence. The distance traveled is significant in this model because it determines the 
locations of ice loading prior to cusp failure over successive icebreaking cycles. The saw¬ 
tooth velocity profile indicates nearly instantaneous momentum change of the icebreaker. 
When all cusps are assumed to fail simultaneously, the model may predict a lower initial 
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velocity limit, below which, the icebreaker is unable to proceed through the ice. Although 
this highlights the motivation for backing and ramming icebreaking operations, where the 
icebreaker cannot proceed continuously, this assumption assumes over predicts the 
propulsive force required to transit the ice for a given initial velocity. To provide a more 
realistic minimum velocity at the start of the icebreaking cycle prior to the necessitation 
of backing and ramming, the velocity profile resulting from the failure of a single ice 
cusp was maintained. This is also shown in Figure 8 as a saw-tooth velocity profile 
having significantly less cyclical amplitude than that predicted with the simultaneous 
cusp failure assumption. 


Vladivostok: Failure of Single v. All Cusps Velocity Profile at 



0.3964 
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Distance Traveled [m] 


- All Cusps 

Single Cusp 


Figure 8. Cyclical cusp failure transient velocity profile. 


2. Influence on Total Resistance 

The steady-state velocity was converged upon for a resistance range of 534 kN, as 
expected for zero-velocity by the Lindqvist equation, to 1400 kN. As shown in Figure 9, 
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the continuous cusp failure assumption exactly overlaid the Lindqvist calculation. The 
cyclical cusp failure assumption, on the other hand, was unable to match neither the 
slope, nor the slope-intercept of the Lindqvist calculation. 


Vladivostok: Continuous v. Cyclical Numerical Ice 
Resistance Calculation 



Velocity [m/s] 


- Lindqvist Calculation 

♦ Continuous 
▲ Cyclical 


Figure 9. Comparison of continuous and cyclical cusp failure with icebreaker 

Vladivostok raw data and Lindqvist resistance predictions versus velocity. 


The cyclical cusp failure calculation error was not expected, because this should 
have perfectly matched the Lindqvist calculation line. Flowever, the average error 
between this numerical calculation to the values expected from the Lindqvist calculation 
over the resistance range was less than three percent, as shown in Figure 10. The likely 
cause of this error is from processing limitations. To assist this argument, the steady state 
velocity was obtained with a MATLAB® code included in Appendix B. Inspection of the 
command window output indicated that computational errors developed from calculating 
the square root term of the velocity iterates. 
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Figure 10. Comparison of continuous and cyclical cusp failure resistance-velocity 
relationship against icebreaker Vladivostok raw data. 


D. CONCEPTUAL TEST DETERMINATION 

The hybrid methodology is most easily evaluated when coupling of the resistance 
components is assumed negligible. By assuming zero velocity, thereby negating the 
rotation and hydrodynamic resistance components, the hybrid method takes a form 
resembling the Lindqvist method. To test the concept of multi-sourced superposition, a 
resistance component was obtained from a second model, and compared with the 
Lindqvist component. The bending resistance component was chosen for comparison 
since the Lindqvist approximation has a non-zero value at zero-velocity and since it is the 
only resistance component that is not localized to the hull-ice contact zone. Therefore, 
testing the resistance component has the three-fold purpose as a proof of concept for the 
hybrid method, cross-examination of the Lindqvist bending resistance, and to provide the 
ability to quantify the effects of damage to the ice analytically. 
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E. 


FINITE ELEMENT MODEL 


Given the zero-velocity assumption to reduce potential analytical errors, the 
bending resistance was obtained from a finite element interpretation of the icebreaker 
providing a static distributed load to the ice sheet at the bow of the icebreaker, as shown 
in Figure 11. 



Figure 11. Illustration of ice loading, cusp length, and transformed cusp length 
determination in static finite-element ice bending model. 

Similar to existing analytical methods, the ice sheet is assumed homogeneous and 
isotropic, with no spatial or time dependent properties, and supported on a uniform elastic 
foundation with stiffness equal to the product of sea water density and gravitational 
acceleration. The finite element model was generated using ANSYS Static Structural 
toolbox. Given the limitations of finite element modeling, it was desirable to select 
dimensions for the ice sheet such that the ice behaved as a semi-infinite plate consistent 
with prior analytical methods and that the ice had a sufficient number of nodes elements 
through the thickness of the ice sheet to eliminate numerical errors as the ice sheet was 
deformed. Assuming the icebreaking process is symmetrical about the centerline of the 
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icebreaker, only one half of the ice sheet was modeled to increase calculation speed. The 
dimensions of this half sheet of ice were arbitrarily selected as 100 m long by 35 m wide. 

Ice was removed from one edge along the length of the ice sheet inwards to 
approximate the area occupied by the half-breadth of the icebreaker in question. In addition, 
this cut assumed a triangular bow fonn, according to the average waterline angle, to roughly 
approximate the influence of this parameter on the bending resistance. The ice thickness was 
varied through a realistic range of 0.1 m to 2.0 m. The ice sheet was bounded on the aft, 
forward, and outboard edges by the fixed boundary condition. The symmetric boundary 
condition was used along the centerline edge of the ice sheet. A uniform elastic support was 
applied on the underside of the ice sheet. A distributed load was applied at the hull contact 
area to roughly estimate the loading on the ice. To improve accuracy of the solution, mesh 
refinement was perfonned in the vicinity of the hull contact area using three-dimensional 
tetrahedral elements. The solution output was the bending failure stress gradient versus the 
distributed load acting on the ice, included as Figure 12. 


ANSYS 

Noncommercial use only 


MM 


Figure 12. FEM static structural bending failure of the ice modeled in ANSYS. 



The distance of ice failure ahead of the icebreaker bow was measured at multiple 
points with respect to the icebreaker breadth. This distance was obtained from multiple 
ice thicknesses and a one order magnitude range for the ice bending stress, from 0.1 MPa 
to 1 MPa. 
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F. 


ANALYSIS 


As both the loading force and distance of broken ice ahead of the icebreaker are 
obtained from the finite element model, this is simply related to the existing analytical 
formulations for the bending resistance. Two obstacles develop, however, since analytical 
ice bending failure models neither represent the force as a distributed load, nor do they 
represent the distance ice is broken ahead of the icebreaker. The latter of the two is 
accounted for since it is the transformed cusp length. The implications of distributed 
loading, on the other hand, required further insight. As expected, without accounting for 
this inconsistency, the relation between the finite element and Lindqvist predictions for 
the bending resistance was not clear across the range of ice thicknesses. Critical 
observation of the Lindqvist method revealed that it assumed both cusp length and width 
to be equal independent of ice thickness. However, other literary sources indicated that 
the length and width were in fact dependent on ice thickness. Since Lindqvist 
approximates the bending resistance as the sum of point loads, spaced apart by the cusp 
width, this total resistance component is proportional to a cusp width dependent 
approximation. Figure 13 illustrates the load distribution and cusp failure geometry 
assumed as a single icebreaking cycle. 
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Figure 13. Illustration of ice loading and cusp failure geometry from real-world 

observation, the Lindqvist 2-D orthogonal ice cusp failure assumption, and 
the static structural finite element ice failure assumption. 


Therefore, the finite element approximation for the average bending resistance 
force is proportional to the Lindqvist value. 



= A 

F v n 



L / J 

FEM 

l 


Lindqvist 


(4.2) 


where A is the proportionality constant. 

The average bending force predictions obtained by finite element method and 
Lindqvist were compared according to parameter variation through a reasonable range of 
values. In this comparison, the Young’s modulus, bending failure stress, and thickness of 
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the ice, as well as the breadth of the vessel, are clearly defined in the Lindqvist equation, 
whereas the influence of these parameters is not clear within the finite element result. The 
range of parameter variation is provided in Table 7. 


Table 7. Parameter variation range 


Parameter 

Min 

Max 

E 

0.2 Gpa 

20 Gpa 

o b 

0.1 Mpa 

1.0 Mpa 

H 

0.1m 

2.0 m 

B 

10 m 

30 m 


Due to the enhanced ability of the finite element model to approximate the cusp 
failure geometry, some error was expected between the finite element and Lindqvist 
average force predictions. The error was linearized in all cases by multiplying the error 
by the variable parameter to some power. This linearization defined the proportionality 
constant, A, between the finite element and Lindqvist average force predictions. The 
average bending force obtained from the finite element result was then directly applied as 
a substitute to the Lindqvist value. 


r m- p + b' 

V 100 p N , 


(4.3) 


Although this form of analysis has limitations, it generally had good agreement for the 
given parameters. The slope and slope-intercept of the linearized relative error was 
obtained and provided in Table 8. 


Table 8. Parameter modification variables for proportionality constant, A. 


P 

N 

m 

b 

E 

0.095 

4.00E-08 

-602.74 

ob 

0.185 

-0.006 

-41.42 

H 

-0.015 

40.406 

-70.215 

B 

0.9 

-69.964 

968.6 


All graphs in this section are presented together in Appendix C as a quick 
reference to this statement. 
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1. Young’s Modulus 

The Young’s Modulus was varied from 0.2 GPa to 20 GPa to accommodate the 
wide range of values presented in literature, since this variable is expected to change 
significantly in the expected operating environment. 

a. Uncorrected Relative Error 

The relative error between the average bending force obtained via the 
finite element and Lindqvist models are provided in Figure 14, showing nonlinear 
dependence on parameter variation. 



Figure 14. Uncorrected relative error of the finite element average bending force v. 

Lindqvist when varying the ice Young’s modulus. 
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b. Linearized Relative Error 

The data was linearized by multiplying the error by the Young’s modulus 
to the power of 0.095, producing a best fit line with an R" value of 0.998, as shown in 
Figure 15. 


Linearized Relative Error; 
Ice Young's Modulus Variation 



-700 

Young's Modulus [Pa] 


y = 4E-08x- 602.74 
R 2 = 0.998 


- FEM Error 

-Linear (FEM Error) 


Figure 15. Linearized relative error of the finite element average bending force v. 

Lindqvist when varying the ice Young’s modulus. 


c. Corrected Average Bending Force 

The best fit error prediction line was then used to correct the finite element 
prediction to the Lindqvist value according to Equation 4.3. The corrected average 
bending force obtained by finite element method is compared with the Lindqvist 
prediction in Figure 16. The data sets are consistent, as expected from this analysis. 
However, the slope-intercept of the finite element corrected value appears to be 
significantly greater than that predicted by Lindqvist. Since the average bending force 
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curves are nearly vertical as the Young’s modulus approaches zero, this error was 
expected by numerical processing limitations. 



Figure 16. Corrected finite element average bending force v. Lindqvist when 

varying the ice Young’s modulus. 


d. Corrected Relative Error 

To finalize the parameter variation of Young’s modulus, the relative error 
of the corrected finite element prediction for the average bending force was compared 
with the Lindqvist value, as shown in Figure 17. As stated before, significant errors were 
expected as the Young’s modulus approached zero. However, if this region is neglected, 
the finite element obtained force nearly matches the Lindqvist prediction with a 
maximum error of 11.3% and a minimum of 1.5%. 
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Corrected Relative Error; 
Young's Modulus Variation 
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Figure 17. Corrected relative error of the finite element average bending force v. 

Lindqvist when varying the ice Young’s modulus. 

2. Bending Failure Stress 

The error between the Lindqvist and finite element average bending force 
predictions was obtained for the bending stress range of 0.1 MPa to 1 MPa to be 
consistent with the range of values obtained during in-situ measurement. 

a. Uncorrected Relative Error 

The relative error between the average bending force obtained via the 
finite element and Lindqvist models are provided in Figure 18, showing nonlinear 
dependence on parameter variation. 
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Figure 18. Uncorrected relative error of the finite element average bending force v. 

Lindqvist when varying the ice bending failure stress. 


b. Linearized Relative Error 

The data was linearized by multiplying the error by the bending failure 
stress to the power of 0.185, producing a best fit line with reasonable agreement to the 
errors presented, as shown in Figure 19. 
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Figure 19. Linearized relative error of the finite element average bending force v. 

Lindqvist when varying the ice bending failure stress. 


c. Corrected Average Bending Force 

The best fit error prediction line was then used to correct the finite element 
prediction to the Lindqvist value according to Equation 4.3. The corrected average 
bending force obtained by finite element method is compared with the Lindqvist 
prediction in Figure 20. The data sets are consistent, as expected from this analysis. 
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Figure 20. Corrected finite element average bending force v. Lindqvist 
when varying the bending ice failure stress. 


d. Corrected Relative Error 

To finalize the parameter variation of the bending failure stress, the 
relative error of the corrected finite element prediction for the average bending force was 
compared with the Lindqvist value, as shown in Figure 21. Although the error curve is 
nonlinear, the maximum and minimum errors are 2.3% and 0.3%, respectively. 
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Corrected Relative Error; 

Ice Bending Failure Stress Variation 
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Figure 21. Corrected relative error of the finite element average bending force v. 

Lindqvist when varying the ice bending failure stress. 

3. Ice Thickness 

The error between the Lindqvist and finite element average bending force 
predictions was obtained for ice thickness ranging from 0.1 m to 2.0 m. This variable is 
expected to change significantly in the expected operating environment. 

a. Uncorrected Relative Error 

The relative error between the average bending force obtained via the 
finite element and Lindqvist models are provided in Figure 22, showing limited if any 
dependence on parameter variation. 
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Uncorrected Relative Error; 
Ice Thickness Variation 
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Figure 22. Uncorrected relative error of the finite element average bending force v. 

Lindqvist when varying the ice thickness. 

b. Linearized Relative Error 

Although the prediction appears to fit well without correction, the data 
was linearized by multiplying the error by the bending failure stress to the power of - 
0.015, producing a best fit line which mitigated residual errors, as shown in Figure 23. 
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Linearized Relative Error; 
Ice Thickness Variation 
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Figure 23. Linearized relative error of the finite element average bending force v. 

Lindqvist when varying the ice thickness. 


c. Corrected Average Bending Force 

The best fit error prediction line was then used to correct the finite element 
prediction to the Lindqvist value according to Equation 4.3. The corrected average 
bending force obtained by finite element method is compared with the Lindqvist 
prediction in Figure 24. The data sets are consistent, as expected from this analysis. 
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Corrected Average Bending Force; 
Ice Thickness Variation 



Figure 24. Corrected finite element average bending force v. Lindqvist when varying 

the ice thic kn ess. 


d. Corrected Relative Error 

To finalize the parameter variation of the ice thickness, the relative error 
of the corrected finite element prediction for the average bending force was compared 
with the Lindqvist value, as shown in Figure 25. Despite the curve presented in the error 
plot, the maximum and minimum errors were 1.0% and 0.3%, respectively. 
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Corrected Relative Error; 
Ice Thickness Variation 



Figure 25. Corrected relative error of the finite element average bending force v. 

Lindqvist when varying the ice thickness. 

4. Vessel Breadth 

The error between the Lindqvist and finite element average bending force 
predictions was obtained for vessel breadth in the range of 10 m to 30 m, consistent with 
most icebreaking vessels. 

a. Uncorrected Relative Error 

The relative error between the average bending force obtained via the 
finite element and Lindqvist models are provided in Figure 26, showing nonlinear 
dependence on parameter variation. 
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Figure 26. Uncorrected relative error of the finite element average bending force v. 

Lindqvist when varying the vessel breadth. 


b. Linearized Relative Error 

The data was linearized by multiplying the error by the bending failure 
stress to the power of 0.9, producing a best fit line with reasonable agreement to the 
errors presented, as shown in Figure 27. 
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Figure 27. Linearized relative error of the finite element average bending force v. 

Lindqvist when varying the vessel breadth. 


c. Corrected Average Bending Force 

The best fit error prediction line was then used to correct the finite element 
prediction to the Lindqvist value according to Equation 4.3. The corrected average 
bending force obtained by finite element method is compared with the Lindqvist 
prediction in Figure 28. The data sets are consistent, as expected from this analysis. 
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Figure 28. Corrected finite element average bending force v. Lindqvist when varying 

the vessel breadth. 


d. Corrected Relative Error 

To finalize the parameter variation of vessel breadth, the relative error of 
the corrected finite element prediction for the average bending force was 
compared with the Lindqvist value, as shown in Figure 29. Despite the curve 
presented in the error plot, the maximum error was 0.11% and the minimum error 
was 0.09%. 
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Corrected Relative Error; 
Vessel Breadth Variation 
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Figure 29. Corrected relative error of the finite element average bending force v. 

Lindqvist when varying the ice vessel breadth 
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V. DAMAGED ICE SHEET RESISTANCE PREDICTION 


A. TEST SPECIFICATION 

Following the development of the hybrid bending resistance calculation, the finite 
element model of the ice sheet was altered to represent specific damage patterns prior to 
interaction with the icebreaker hull. This modified ice sheet was then loaded with the 
icebreaker hull, and this load was substituted for the analytical bending resistance 
component in the Lindqvist calculation. As the substitution errors are relatively small for 
the range of values selected in this initial study, the variation in total ice resistance 
resulting from the damage to the ice sheet contains minimal error within the limitations of 
the static structural finite element model. Although there are infinite ice damaging 
strategies that may be implemented, this study only entertains the concept of linear cuts 
travelling parallel to the icebreaker at three locations along the breadth of the vessel; at 
the vessel centerline, at lm outboard of the extreme breadth, and 2m outboard of the 
extreme breadth. The cuts were defined as either partial penetration through 50% of the 
ice sheet, or as full penetration cuts and were assumed to occur at some distance well 
ahead of the icebreaker. In total, 8 specific damage patterns were tested. The tests are 
identified in Table 9. 


Table 9. Ice damage test case identification 


Test Case 

Location 

Depth 

1 

Side [lm] 

50% 

2 

Side [lm] 

100% 

3 

Side [2m] 

50% 

4 

Side [2m] 

100% 

5 

Center 

50% 

6 

Center 

100% 

7 

Side [lm] and Center 

50% 

8 

Side [lm] and Center 

100% 
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B. 


RESULTS 


The ice was damaged according to the patterns defined in the test specification, 
and this was accounted for as an averaged bending force. The finite element bending 
failure stress contours are included in this section to illustrate the location of damage and 
the expected failure regions. The obtained bending resistance component for damaged ice 
relative to undamaged ice is included in Table 10. 


Table 10. Relative bending resistance of undamaged ice versus damaged ice 

by test case. 


Test Case 

Location 

Depth 

Rb [%] 

1 

Side [lm] 

50% 

132.40 

2 

Side [lm] 

100% 

11.16 

3 

Side [2m] 

50% 

119.54 

4 

Side [2m] 

100% 

147.52 

5 

Center 

50% 

159.23 

6 

Center 

100% 

24.63 

7 

Side [lm] and Center 

50% 

68.17 

8 

Side [lm] and Center 

100% 

6.00 


1. Side Cut [lm] 

This cut involves a cut lm outboard of the extreme breadth of the icebreaker, 
running parallel to the vessel’s longitudinal dimension in the finite element model. 

a. Partial 

The cut was initiated at the upper surface of the ice through 50% of its 
thickness. The ice directly ahead of the icebreaker failed similarly to the undamaged ice 
model; however, the stress contour is visible at the cut location. The finite element 
solution to this damage pattern for the bending failure stress contour is included as Figure 
30. 
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Figure 30. Bending failure stress contour for a partial side cut located lm outboard 

of the extreme breadth of the icebreaker. 


b. Full 

The cut was initiated at the upper surface of the ice through the entire 
thickness, constituting the removal of the entire thickness of the ice at the location of the 
cut. The ice failure region differs considerably from the undamaged ice model, indicating 
stress concentrations at both the forward most extent of the cut and at the icebreaker 
shoulder. The finite element solution to this damage pattern for the bending failure stress 
contour is included as Figure 31. 






».00«(m) 


Figure 31. Bending failure stress contour for a full penetration side cut located lm 

outboard of the extreme breadth of the icebreaker. 
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2. Side Cut [2m] 

This cut involves a cut 2m outboard of the extreme breadth of the icebreaker, 
running parallel to the vessel’s longitudinal dimension in the finite element model. 


a. Partial 

The cut was initiated at the upper surface of the ice through 50% of its 
thic kn ess. The ice directly ahead of the icebreaker failed similarly to the undamaged ice 
model, however, the stress contour is visible at the cut location. The finite element 
solution to this damage pattern for the bending failure stress contour is included as Figure 
32. 
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Figure 32. Bending failure stress contour for a partial side cut located 2m outboard 

of the extreme breadth of the icebreaker. 
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b. Full 

The cut was initiated at the upper surface of the ice through the entire 
thickness, constituting the removal of the entire thickness of the ice at the location of the 
cut. The ice failure region differs considerably from the undamaged ice model, indicating 
stress concentrations at both the forward most extent of the cut and at the icebreaker 
shoulder. The finite element solution to this damage pattern for the bending failure stress 
contour is included as Figure 33. 
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Figure 33. Bending failure stress contour for a full penetration side cut located 2m 

outboard of the extreme breadth of the icebreaker. 


3. Center Cut 

This cut involves a cut at the centerline of the icebreaker, and running parallel to 
the vessel’s longitudinal dimension in the finite element model. 

a. Partial 

The cut was initiated at the upper surface of the ice through 50% of its 
thic kn ess. The ice directly ahead of the icebreaker failed similarly to the undamaged ice 
model, however, a high stress region extends forward of the stem of the icebreaker at the 
cut location. The finite element solution to this damage pattern for the bending failure 
stress contour is included as Figure 34. 
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Figure 34. Bending failure stress contour for a partial penetration centerline cut located 

at the stem of the icebreaker. 


b. Full 

The cut was initiated at the upper surface of the ice through the entire 
thickness, constituting the removal of the entire thickness of the ice at the location of the 
cut. The ice directly ahead of the icebreaker failed similarly to the undamaged ice model; 
however stress concentrations occurred at both the forward most extent of the cut and at 
the icebreaker shoulder. The finite element solution to this damage pattern for the 
bending failure stress contour is included as Figure 35. 



Figure 35. Bending failure stress contour for a full penetration centerline cut located 

at the stem of the icebreaker. 
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4. 


Side and Center Cut 


This cut involves both a cut lm outboard of the extreme breadth of the icebreaker 
and at its centerline, running parallel to the vessel’s longitudinal dimension in the finite 
element model. 

a. Partial 

The cut was initiated at the upper surface of the ice through 50% of its 
thickness. The ice directly ahead of the icebreaker failed similarly to the undamaged ice 
model; however, high stress regions develop at the shoulders and directly forward of the 
stem of the icebreaker at the cut locations. The finite element solution to this damage 
pattern for the bending failure stress contour is included as Figure 36. 



Figure 36. Bending failure stress contour for partial penetration side and centerline 

cuts, located at the extreme breadth and stem of the icebreaker, respectively. 


b. Full 

The cut was initiated at the upper surface of the ice through the entire 
thickness, constituting the removal of the entire thickness of the ice at the location of the 
cuts. The ice failure region differs considerably from the undamaged ice model; 
indicating stress concentrations at both the forward most extent of the cut and at the 
icebreaker shoulder. The finite element solution to this damage pattern for the bending 
failure stress contour is included as Figure 37. 
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Figure 37. Bending failure stress contour for full penetration side and centerline cuts, 
located at the extreme breadth and stem of the icebreaker, respectively. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

1. Hybrid Method 

A hybrid method was developed such that ice resistance data from multiple 
sources may be cross-referenced or combined under common framework to better 
describe icebreaking in the analytical context. Within this framework, a numerical 
approximation for the bending ice resistance component was obtained via finite element 
model and related to existing analytical model with empirical basis in order to better 
understand the underlying physics as well as to improve the accuracy of existing ice 
resistance models. Upon achieving low error for given predefined influential parameters, 
the finite element ice sheets were damaged according to specified damage patterns, and 
this damage was accounted for as a decrease in the total ice resistance by using the 
Lindqvist analytical model. 

2. Damaged Ice Sheet Resistance 

The finite element model of the ice sheet was altered according to eight damage 
scenarios to relate any change in the bending stress failure contour to the bending 
resistance component. The method presented in this study indicates that the bending 
resistance is significantly reduced when the ice is damaged by either full penetration cuts, 
or cuts extending ahead of the icebreaker stem. The reduction of bending resistance for 
the centerline cut is particularly interesting because this is also the location of the 
breadth-independent stem crushing component. It may be possible to reduce the stem 
crushing component as well as the bending resistance for centerline cuts in the ice. Since 
the stem crushing and bending resistance components are likely coupled, this study was 
unable to account for any reduction in the stem crushing resistance while solving for the 
bending resistance. The present finite element static structural model could not account 
for stress gradient or stress rate failure criteria for the ice sheet because these values are 
not well represented in literature. 


63 



B. FUTURE WORK 

To proceed with this work, both the time and size dependent failure properties of 
sea ice must be obtained. Without this information, it was impossible to numerically 
model the transient failure of ice and test against empirical data. In addition, further tests 
must be perfonned to quantify resistance component coupling in the hybrid model 
framework so that analytical methods may be further refined. If, for instance, the 
coupling between the stem crushing and bending resistance components are resolved, the 
effects of ice damage on the stem crushing component may be pursued. 
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APPENDIX A. LINDQVIST SENSITIVITY ANALYSIS 


Each parameter in the Lindqvidst resistance formula was varied across a 
reasonable range of values, with respect to the vessels used by Lindqvist to formulate the 
method. The figures below illustrate the effect of this variation on the individual 
resistance components. 



65 























Draft Sensitivity 

600 

500 

1, 400 

v 

c 300 — ■'—' 

re 

•R 200 _ Submergence 

oc 

100 

0 -I - 1 - 1 - 1 - 1 

0 5 10 15 20 

Draft [m] 


Waterline Entrance Angle (a) 
Sensitivity 



0 10 20 30 40 50 60 70 80 90 


Waterline Entrance Angle (a) [degrees] 


♦ Stem Crushing 

- -Bending 

Submergence 


66 



















V 

cc 


8000 

6000 

4000 

2000 

0 

-2000 

-4000 

-6000 

-8000 


Stem Angle Sensitivity 


<> — ♦ ♦ t ♦ 


T-1-1-1-1 


0 10 20 30 40 50 60 70 80 90 


Stem Angle (y) [degrees] 


♦ Stem Crushing 


67 























Stem Angle Sensitivity 



Stem Angle (y) [degrees] 


- —Bending 


Stem Angle Sensitivity 



Stem Angle (y) [degrees] 


Submergence 


68 






















Ice Thickness Sensitivity 



•• Stem Crushing 
— Bending 
— Submergence 


Bending Failure Stress Sensitivity 



Bending Failure Stress [kPa] 


. Stem Crushing 

- -Bending 


69 

























Young's Modulus Sensitivity 


160 
140 
^ 120 
— 100 


<u 

C£ 


80 

60 

40 

20 

0 


-V 


■v 




5E+09 


1E+10 


Young's Modulus [Pa] 


- —Bending 


70 



























Poisson's Coefficient Sensitivity 



- -Bending 


Poisson's Coefficient 


Sea Water/Ice Density Differential 
Sensitivity 



Submergence 


71 





















THIS PAGE INTENTIONALLY LEFT BLANK 


72 



APPENDIX B. NUMERICAL ICEBREAKING VELOCITY 

MATLAB CODE 


close all 
clear all 
clc 

% Icebreaker Vladivostok 
% Parameters 

L = 112; 

B = 2 3.5; 

T = 9.5; 

alpha = 18; % waterline entrance angle 
gamma =24; % stem angle 

zeta = 56; % effective ice contact angle 
mu = 0.16; 

H = 0.58; 

Hsnow = 0.28; 

E = 2E9; 
nu = 0.3; 
sigmab = 500; 
sigmac = 3430; 
rhow = 1026; 
deltarho = 109; 

m = 5000; 
g = 9.81; 

Fpmin = 0; 

Fpmax = 1400; 

Fv = 84.1; 

Rc = 67.25344; 

Rb = 105.383; 

Rs = 360.9507; 

lc = (E*H A 3/(12*(l-nu A 2)*rhow*g)) A 0.25; 

1 = lc/3; 

lx = 1/sind(alpha); 
lxn = 10; 
lxi = lx/lxn; 

a = Fv/(sigmac*10 A 3*B); 
ax = a/sind(alpha); 
axn = 10; 
axi = ax/axn; 

s = ax+lx; 
sn = axn+lxn; 

n = B/(l*sind(alpha)); 
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Fbmax = Rb*(ax+lx)/(2*ax); 


for k = Fpmin:Fpmax 
Fp = k; 

Fprop(k-Fpmin+1,1) = Fp; 

i = 1; 

Vn = 1.5; 

Vs = 0; 

Vsf = 1; 

while Vsf-Vs > 0.05 
for j = 1:axn 
Vs = Vn; 

Fb = Fbmax/axn*j; 

Rci = Rc*(1+1.4*Vn/sqrt(g*H)); 

Rbi = Fb*(1+1.4*Vn/sqrt(g*H)); 

Rsi = Rs*(1+9.4*Vn/sqrt(g*L)); 

Vi (j,i) = Vn; 

Vo = Vi(j,i); 

Vx = 2*axi/m*(Fp-Rci-Rbi-Rsi)+Vo A 2; 
Vi(j+l,i) = sqrt(Vx); 

Vn = Vi(j +1,i); 

V(j+(i-1)*(sn+1)) = Vn; 

end 

for j = axn+l:sn+l 
Vi (j,i) = Vn; 

Fb = 0; 

Rci = Rc*(1 + 1.4*Vn/sqrt(g*H)) ; 

Rbi = Fb*(1+1.4*Vn/sqrt(g*H)); 

Rsi = Rs*(1 + 9.4*Vn/sqrt(g*L) ) ; 

Vo = Vi(j,i); 

Vx = 2*lxi/m*(Fp-Rci-Rbi-Rsi)+Vo A 2; 
Vi(j+l,i) = sqrt(Vx); 

Vn = Vi(j +1,i); 

V(j+(i-1)*(sn+1)) = Vn; 

end 

Vsf = Vn; 
i = i +1; 

end 

Vss(k-Fpmin+1,1) = k; 

Vss(k-Fpmin+1,2) = mean(Vi(:,i-1)); 

end 

plot(V) 

%Lindqvist = [0 534;4.251 1400]; 

%plot(Vss(:,2),Vss(:,1)) 

%hold on 

%plot(Lindqvist(:,1),Lindqvist(: , 2) ) 
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APPENDIX C. HYBRID METHOD DEVELOPMENT ERROR 


Linearized Relative Error 
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